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1. Introduction 

Gamma-ray bursts are the most powerful explosions in the Universe. (For reviews see 
Piran 1999 and Meszaros 2002.) A crucial advance in understanding gamma-ray bursts 
began with the discovery of “afterglows” , starting with BeppoSax observations in the soft 
X-ray band of GRB 970228 (Costa et al. 1997; Wijers, Rees, k Meszaros 1997). (For a 
review of GRB afterglows see van Paradijs, Kouveliotou, & Wijers 2000.) If GRBs were 
isotropic, then the measured redshifts would imply total explosion energies of ~ 10 52 — 10 54 
ergs (Frail et al. 2001). Theoretical work on relativistic jet expansion, however, shows 
that one expects a steepening in the decay light curve if one is looking down the axis 
of a jet as the flow decelerates from a bulk Lorentz factor < 9 to 7 -1 > 0, where 9 
is the jet beaming angle (e.g., Rhoads 1997, Sari, Piran, k Halpern 1999, Panaitescu & 
Meszaros 1999, Panaitescu & Kumar 2001ab, 2002). One does in fact see such steepenings 
in the light curves (e.g., Stanek et al. 1999; Harrison et al. 1999). Prior to the time when 
7 -1 ~ 9 the expansion is effectively “spherical” from the observer’s viewpoint because the 
relativistic beaming is narrower than the jet itself. In other words, if the GRB emission 
were coming from one spot on a large, relativistically expanding sphere, aimed directly at 
the observer, the observer would not see any emission from the other parts of the sphere. 
After the time when 7 -1 ~ 6 the observer can “see” the entire jet, and a faster rate of 
decline in the luminosity is predicted. A separate issue that we will address in this work is 
the sideways or lateral expansion of the jet as the increasing solid angle of the jet enables 
a larger fraction of the circumstellar medium (CSM) surrounding the progenitor star to 
be intercepted and provide decelerating gas. Several groups have claimed that this leads 
to a faster (exponential) decrease in 7 , which acts as an additional agent to diminish the 
amplitude of the relativistic beaming. 

The concept of a “break” in the afterglow light curves occurring when 7 -1 ~ # has 
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the jet, and (iii) a deceleration phase during which the accumulated mass forces a rapid 
decrease in 7. Previous studies using analytical methods divided the evolution into regimes 
defined by some ordering of distance, velocity, or energy scales (see Table 3 of Piran 1999 
and Figure 7 of Meszaros 2002). The main findings that seem to be common to all studies 
are that the break in the decaying light curve d log L(t)/d log t occurs roughly when the 
deceleration has decreased 7 to roughly the reciprocal of the jet beaming angle, and that 
the subsequent decrease in 7 is roughly exponential with distance. Also, Panaitescu & 
Meszaros (1999) present axisymmetric calculations to study the combined effects of the 
transition from 7 -1 <9 to 7 -1 > 9 and the lateral jet expansion. 

In other subdisciplines of astronomy the use of relativistic hydrodynamics codes has 
been standard for some time. For example, workers studying extragalactic jets have used 
such codes to continuously inject a collimated supersonic beam into a surrounding medium, 
usually under the assumption of pressure equilibrium, (see, e.g., Norman et al. 1982 for 
a thorough discussion). In the context of GRBs, work has been done using 2D and 3D 
relativistic hydro codes to consider the evolution of the GRB jet as it propagates through 
the envelope of the progenitor star, up to the point where it breaks out of the stellar surface 
and produces the prompt GRB emission (Zhang, Woosley, & MacFadyen 2003, Zhang, 
Woosley, & Heger 2003). In this work we consider the evolution covering the afterglow time 
(i.e, after the period considered by Zhang et al. 2003ab). We utilize a three dimensional 
relativistic hydrodynamical code to study the propagation of an initially ultrarelativistic 
blob into a dense CSM. We study the spatial spreading of the blob both along the direction 
of propagation and orthogonal to it, as well as the evolution of 7 in space and time. We also 
calculate afterglow light curves, taking a simple prescription in which the local emissivity 
scales with the local pressure. 
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DB02 present a suite of test results for 1, 2, and 3 dimensions. We have reproduced 
these tests, and now show the results of a 3D spherical expansion that results in higher 7 
values than the tests discussed in DB02. We utilize a cubical grid of 100 x 100 x 100 points, 
with each side having length unity. 2 Initial conditions are that everywhere = 0, p = 2, 
and p = 3 x 10 5 within a radius 0.1 of the center of the cube, and p = 1 outside. This 
is the 3D analog of the ID “piston” problem. (An identical test with less extreme initial 
conditions is shown and discussed by Hughes, Miller, Sc Duncan 2002, see their Fig. 3.) 
The enormous over-pressure launches a relativistic blast wave at t = 0. Figure 1 shows the 
evolution of p , v x , p, and 7 for a slice taken along the x— axis. The density inside the sphere 
becomes small for late times. To avoid numerical instabilities in this region we utilize the 
smoother LLF fluxes for small densities (p < 0.0225), and the HLL fluxes elsewhere. We 
show the initial configuration plus 10 time slices from a run with 360 total time steps, 
taking a Courant number of 0.05. After the run begins one can see the development of a 
strong spherical outflow from the center of the cube. The Lorentz 7 factor has increased to 
~ 20 — 30 in the outer parts of the expanding shell by the end of the run. The run was 
halted before the expansion reached the edge of the grid. For the last three times steps, a 
kink develops in the v x profile at the center of the sphere because of the prolonged spherical 
outflow from localized point. Figure 2 shows the corresponding evolution for pure LLF 
fluxes. All variables now show a smoother evolution and 7 — *• 20 at large radii for late times 
— somewhat less than for the HLL fluxes. Conservation of rest-mass energy M = Sj(pj7j) 
and energy E = 'E i (e i 'yf + #(7? — 1)) (where the specific energy e< = p x + Pi/(T — 1)) 

2 DB02 quote a CPU run time of ~ 2 minutes per Runge-Kutta integration for a 3D 

problem with 100 3 * nodes, utilizing the third order reconstruction for the L and R values of 

the primitives, and using a 1 GHz Linux PC; our run time of ~ 70 s is consistent with the 

1.7 GHz Linux PC used in this work. 



-9- 


relaxes to one in which Pbiob/pcsM — 10, therefore the effective density contrast for the early 
run is ~ 200. This is less than expected astrophysically, but required in our computations 
in order to see significant deceleration of the blob by the time it reaches the end of the grid 
x = 5. In other words, we must telescope the evolution into the finite dimensions of our 
grid. Prom conservation of momentum, the condition for significant deceleration is that the 
total swept-up mass-energy roughly equals that in the initial blob. 

Experimentation using HLL fluxes shows that strong internal shocks almost 
immediately create large-amplitude sawtooth p variations within the blob, leading to noisy 
results. Therefore we utilize the smoother LLF fluxes in this work. We follow the evolution 
of the blob in terms of both its motion in x and its spreading in y and z. Taking a Courant 
number of 0.25 necessitates ~ 2000 time steps for the blob to reach the end of the grid, 
traveling at v x ~ 1. This is a simple consequence of the box length and grid spacing 
Ax = 0.01, from which it follows that the time to traverse the box is 2000 x 0.25 x 0.01 = 5 
units traveling at v ~ 1. A density contrast of ~ 10 2 is sufficient to see the desired 
deceleration from y" 1 < 9 to 7 -1 > 9 during a run. 

In order to avoid undue complexity in these experiments, that currently are purely 
hydrodynamic and do not yet contain proper prescriptions for emission from bremsstrahlung 
and synchrotron processes, we calculate a simple measure of the emissivity by taking the 
local emission to scale as p, which would be expected roughly for optically thin synchrotron 
emission characteristic of frequencies significantly above the self-absorption frequency. 

We avoid the issues of synchrotron self-absorption (Granot, Piran, & Sari 1999b, 2000; 
Panaitescu & Kumar 2000; Sari & Esin 2001), and of whether the evolution is adiabatic or 
radiative (Panaitescu k, Meszaros 1998b; Sari, Piran, & Narayan 1998) — our calculations 
are adiabatic. We calculate effective “light curves” for observers at various viewing angles 
between 0° and 15° from the center of the jet. The amplification of the photon energy flux is 
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behind this shocked region. Although the initial pressure is small, after the bow shock 
is fully developed we have p S hock/Pshock — 0-3 throughout the subsequent evolution. The 
evolution proceeds in a roughly self-similar manner and the forward/reverse shock system 
is stable, in accord with analytic estimates (Wang, Loeb, & Waxman 2002). 

Figures 4 and 5 show a time series of contour plots in 7 and p that follow the blob 
propagation and radial spreading. The full 500 grid points along the x direction are shown, 
and each frame represents 240 time steps (i.e., A t = 0.6). The leading contours show the 
density enhancement associated with the shock. The p contours indicate the values 0.2, 
0.3, 0.4, and 0.5, while the (trailing) 7 contours, indicate values of 1.25, 2.5, 5, 10, 15, 
and 20. The contours are formed by taking a cut at z — 0 through the x — y plane. The 
initial velocity vector of the blob points toward the top of the plot, with 7 = 25. One sees 
a strong bow shock associated with the maximum in pressure at the point where the blob 
encounters the GSM. The lack of a bow shock in the 7 contours suggests that little or no 
material in advance of the blob is accelerated to a significant bulk Lorentz factor. As the 
evolution progresses the deceleration causes the higher 7 contours to disappear, and those 
that remain become increasingly distorted. The bilateral symmetry of the contours evident 
in Figures 4 and* 5 attests to the power of the third order differencing scheme given in 
DB02, insofar as the 3D model has no enforced symmetry. 

Figure 6 shows the evolution of the total rest-mass energy M (dotted) and total energy 
E (dashed) within the grid, over the course of the conical run. The values of M and E 
have been normalized to their initial values. The curves that extend significantly below an 
abscissa value of unity indicate the values summed over the grid. At the six faces of the 
computational box we continually reset the values of all variables to their initial values, and 
keep track of the differences between those values and the initial ones. As the run progresses 
and more high velocity gas reaches the edge of the computational domain and is extracted, 
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different 7 cuts follow the values 7 cu t = 1.001, 1.01, 1.1, 1.5, 2, 4, 6, and 8. The higher 7 
cuts 7 cut > 2 reveal that a negligible fraction of CSM matter gets accelerated to significant 
values. This is because the curved bow shock shunts material laterally in front of the 
advancing blob, rather than accelerating it up to a significant fraction of the blob’s bulk 
Lorentz factor ~ 10 — 20. The rest-mass energy curves are relatively constant in time up to 
t ~ 2. At roughly this point the total accumulated CSM gas becomes comparable to that 
in the initial blob. If the blob did not spread laterally, this time would occur when the blob 
had plowed through a density ~ 7<nitiai(Pbiob/pcsM), where initial = 20 and Pbiob/pcsM — 10. 
This means that the blob would have shed of order its initial momentum by the time it had 
traveled ~ 200 times its initial length ~ 6x = 0.1, or ~ 20 units. In practice, the initially 
imposed tangential velocity vt leads to the lateral expansion that effectively increases 
the cross section for interaction of the blob. The negligible decrease in < ^v x > for the 
t> 7 >( 0 ) = 0 run demonstrates the importance of this effect. 

Figure 8 shows the lateral expansion of the blob for the conical run. The blob edge is 
computed in each time step by first finding the local maximum in either p , 7, or pj along 
the x— axis, and then stepping laterally to the position at which the background (CSM) 
value of the relevant quantity has been increased by 10% due to expansion of the blob. The 
transient lateral expansion v edge ~ 0.3 — 0.4c lasting until x ~ 1 is unphysical insofar as 
it occurs during the period of adjustment to the initial density and velocity profiles. The 
value of the later spreading rate Vedge ~ 0.1c is basically dictated by the initial v? value 
given to the gas. 

Figure 9 shows the evolution of < -yv x > and 8 as the blob propagates. The solid curve 
in the top panel indicates the same weighted value of < jv x > as shown in Fig. 7. After 
the transient physical conditions associated with the initial state have vanished, one sees a 
period of deceleration associated with the increasing drag force. The dashed curves indicate 
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4. Discussion 

Utilizing 3D relativistic hydrodynamical calculations, we have examined the evolution 
of an expanding relativistic blob of gas intended to be representative of a jet associated with 
ejecta from an extremely energetic event such as a hypernova, that produces a gamma- ray 
burst (Aloy et al. 2000; Tan, Matzner, & McKee 2001; MacFadyen, Woosley, & Heger 
2001, Zhang, Woosley, & Heger 2003, Zhang, Woosley, & MacFadyen 2003). Since these 
are the first such calculations applied to the blob during the time in which the afterglow 
radiation is produced, we have purposely kept them simple in an effort to concentrate on 
the most fundamental aspects of the physics. We restrict our attention to the transition 
from spherical to jetlike expansion that occurs during the time that the Lorentz factor 
becomes less than the reciprocal of the jet spreading angle. 

We have not yet attached specific numbers to our results. From the SRHD equations, 
one sees that the relevant quantities are the ratios of pressure to density, and of distance to 
time. If we specify either one of these two sets of numbers, the other one is also determined. 
The column giving the observed afterglow break time tj in Table 1 of Frail et al. (2001) 
indicates tj ~ 2 d as being representative. For an observer directly on the velocity vector of 
the blob, the time T between the GRB and afterglow 
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where the dominant contribution to the integral comes from later times. Thus the light 
travel time of 1 day is multiplied by ~ 27 2 fterglow ~ 2 x 10 2 , assuming the spherical-to-jetlike 
transition giving the break in the afterglow light curve happens at 7 ~ 10. For the conical 
run, the break in the light curve occurs at at x ~ 2.5. If we designate this point as 
corresponding to a time 2 d, then x = 2.5 translates to ctj ~ 5.2 x 10 15 (200) cm ~ 10 18 
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precursor and continues to propagate, the medium through which it travels should be 
dominated by the density profile left from the precursor’s stellar wind (Chevalier & Li 
1999, 2000, Li & Chevalier 2001). If the density varies as r~ 2 away from the star, the 
mass loss of the wind Mwind — 10 ~ 5 M© yr -1 , and the wind velocity v mnd ~ 10 3 km s _1 , 
the circumstellar density at the point where the afterglow radiation is emitted ~ 10 18 cm 
will be n C sM ~ Mwind/(47rr 2 v wind m p ) ~ 0.3 cm -3 . In the calculation of Zhang et al (2003), 
the density inside the jet is ~ 10 17 cm -3 at ~ 10 12 cm, and p ~ p. If the spreading angle 
of the jet remains roughly constant from ~ 10 12 cm to ~ 10 18 cm (Lithwick & Sari 2001), 
then the density inside the jet at ~ 10 18 cm should be lower by ~ (r 2 /ri) 2 ~ 10 12 , or about 
~ 10 5 cm" 3 . At this point p inside the jet should be negligible compared to p. The density 
contrast between the blob and CSM is greater than what we have assumed in the results 
shown previously. A separate run taking an initial density contrast of 10 6 shows the same 
basic effects as the previous runs, however, namely a lateral expansion rate of ~ 0.1c and 
the non-accretion of CSM gas. 

Our computations lie within the “deceleration phase” discussed by Kobayashi et 
al. (1999, see also Kobayashi &; Sari 2000, 2001). We find a change in the form of the 
luminosity decrease corresponding to the transition between spherical and jetlike expansion. 
The determination of the average spreading angle <6 > is nontrivial because it depends on 
how one does the averaging, and how much of the diffuse, sideways-expanding jet material is 
included in the computation. In Figure 9 we presented cuts for gas possessing 7 > 7 cut = 4, 
5, and 6 as representative of material in the flow that partakes most strongly in producing 
the observed radiation. We do not see a dramatic increase in < 9 > during the deceleration 
phase; the < 0 > value basically reflects the ballistic motion of material following its initial 
vt value. Also, because the deceleration in our problem is forced not by the accumulation 
of gas from the CSM but rather the drag force of the CSM on the blob, we do not see an 
exponential decrease in 7 with distance during the deceleration phase, but rather a decrease 
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During the latter stage, the swept-up mass increases exponentially in time. Panaitescu & 
Meszaros (1999) calculate light curves for observers at varying angles from the jet axis, and 
calculate separately the effects of including and excluding the lateral jet expansion. They 
find that the maxima in the light curves occur substantially later in runs which do not 
take into account the jet broadening (see their Fig. 4). In the analytical models of both 
Rhoads (1999) and Panaitescu & Meszaros (1999), the physics of mass accumulation from 
the CSM is an integral component of the formalism; all mass within the solid angle of the 
expanding shell is assumed to accrete. Workers have applied the results of Rhoads (1999) 
and Sari, Piran, & Halpern (1999) to the afterglow evolution, however the results of Zhang 
et al. (2003) cast doubt on the validity of this exercise, because at the time corresponding 
to the afterglow emission one anticipates that p « p and the lateral spreading rate would 
not be governed by the internal sound speed but rather the ballistic motions of the ejecta 
comprising the blob as they leave the vicinity of the progenitor star. The blob has a large 
internal thermal energy p/p — 10 as it emerges from the progenitor star, and even by the 
time the expanding ejecta have become optically thin to their own emission from internal 
shocks (producing the GRB), one still expects p/p ~ 1. By the (much later) time of the 
afterglow emission, however, the blob would have cooled to the point that p/ p << 1, 
which provides the impetus for our initial condition Pbiob/Pbiob = 10~ 4 . Following the 
arguments of Rhoads (1999) and Sari, Piran, & Halpern (1999), if the lateral spreading 
rate were mandated by the internal sound speed, then in our calculations it should be 
~ \fpJ~P = 0.01c, whereas we find it to be ~ 10 times larger. The spreading rate basically 
reflects our initial vt value. Thus the physics of the lateral expansion is different than in 
Rhoads (1999) and Sari, Piran, & Halpern (1999). How might this result be influenced by 
systematic effects present in our model? One obvious potential shortcoming is the absence 
of cooling. In this work we have assumed an adiabatic gas, whereas in reality one might 
envision the presence of cooling within the shock. This might then reduce the ability of 



- 21 - 


where Uo is the initial value of U and a = (K/3)pcsMC r /M>iob- From the fitting to < U > 
presented in Figure 7, we infer that K ~ 2.4. The specific numerical value for K is probably 
influenced by our numerical resolution. A comparison of the curves labelled “C” and “f(x) n 
in Fig. 7 shows that the functional decrease in < U > with x is reasonably described by a 
cubic, as expected if the cross sectional area increases quadratically with x , or equivalently 
t. Note that we assumed v x = 1 in this exercise, which is a good approximation for the 
evolution of interest. 


5. Conclusion 

The calculations we present are the first 3D relativistic hydrodynamical calculations 
of GRB jet evolution pertinent to the afterglow phase that do not enforce any special 
symmetry (e.g., spherical or axial). We find that (i) the CSM gas does not accrete onto 
the advancing blob, but rather is shunted aside by the bow shock, (ii) the decay light curve 
steepens roughly when one first “sees” the edge of the jet 7 -1 « 9, with this effect being 
strongest for “face-on” observers (confirming previous studies), and (iii) the rate of decrease 
of the x — component of momentum < jv x > is well-characterized by a simple model in 
which the cross sectional area of the blob increases quadratically with laboratory time 
(or distance). The primary impetus for the built-in assumption of accretion of matter in 
previous studies was the influential work of BM76 in which spherical relativistic expansion 
was considered. Accretion of gas onto the relativistically expanding shell is obviously 
justified for spherical expansion, but subsequent GRB workers applied the results to the 
case of the GRB jet, in which a thin wedge of material propagates through a low density 
medium. In such a situation the natural tendency of material in front of the jet is to 
be pushed aside and to form a “channel flow” around the jet, rather than to accrete. A 
separate issue is that workers used the results of Rhoads (1999) and Sari et al. (1999) that 
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<t>, averaged along the axis of the jet, seems questionable. In addition, one would still be 
missing important physical effects, such as the non-accretion of the forward CSM material. 

Two obvious refinements, currently being carried out, are to (i) treat the problem on 
a Lagrangian grid in which the mesh points follow the blob and are adaptively inserted 
in regions with strong gradients, so as to be able to explore regimes in which the density 
contrast between the blob and CSM is much larger, and (ii) include provisions for realistic 
bremsstrahlung and synchrotron physics, in order to produce light curves that can be 
compared directly with observations so as to test different aspects of the theory and thereby 
constrain the allowed parameter space. 

Our sincerest thanks go to Luca Del Zanna who provided key insights into the workings 
of DB02. We also thank John Baker, David Band, Tom Cline, Chris Fragile, Chris Fryer, 
Markos Georganopoulos, Peter Goldreich, Demos Kazanas, Pawan Kumar, Zhi-Yun Li, 
Andrew MacFadyen, Peter Meszaros, Ewald Mueller, Jay Norris, Alin Panaitescu, Sterl 
Phinney, Steve Reynolds, James Rhoads, Steve Ruden, Sabrina Savage, Craig Wheeler, and 
Weiqun Zhang for helpful comments. 
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FIGURE CAPTIONS 

Figure 1. The evolution of a spherical relativistic expansion due to a large over-pressure 
inside a sphere of radius 0.1, adopting HLL fluxes. The evolution encompasses 360 time 
steps, taking a Courant number of 0.05. Shown are the initial conditions plus ten equally- 
spaced time steps taken from a slice along the x — axis depicting the evolution of (i) pressure 
p ( top panel), (ii) density p ( second panel), (iii) v x /c ( third panel), and (iv) Lorentz factor 7 
( bottom panel). The small numbers beside each curve indicate the evolution. One sees an 
expansion of the central over-pressurized sphere into the surrounding medium. The strong 
relativistic outflow peaks at 7 ~ 20 — 30 near the end of the evolution. The general trend 
in which 7 oc r within the expansion can be derived from fundamental principles and is well 
known from classical solutions. At late times the strong decrease in p at the center of the 
sphere, which is a singularity in this test, begins to cause numerical instability. We utilize 
the LLF fluxes for p < 0.0225 and the HLL fluxes for p > 0.0225. 

Figure 2. The evolution of a spherical relativistic expansion due to a large over-pressure 
inside a sphere of radius 0.1, adopting pure LLF fluxes. Quantities shown are the same 
as in Fig. 1. The evolution is nearly identical to that in Fig. 1, except the variations in 
the physical variables are smoother, and 7 ~ 20 at larger radii by the end of the run — 
somewhat smaller than in Fig. 1. 

Figure 3. The evolution of a relativistic blob launched as a small cone traveling along 
the -f-x-axis, initially confined to between x = 0.1375 and x = 0.2 and a maximum radius 
0.025. The small numbers beside each curve indicate the time step. The panels are the 
same as shown in Fig. 1. Snapshots represent conditions along a slice through the center 
of the rectangular grid, taken every 120 time steps (At = 0.2). For this trial and all that 
follow we utilize a Courant number of 0.25 and pure LFF fluxes. The initial bulb Lorentz 
factor 7 was set to 25, and the initial full width flaring angle to 0.07 radian as measured by 
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to the conical run given by 21 — 0.15(z + 0.45) 3 . Although the rest mass energy curves for 
Tcut = 1001 and 1.01 show a sharp increase (up until t ~ 1.7 when material starts to leave 
the edges of the grid), the near constancy of the 7 cut = 2 curve and the slight decline for 
higher 7 cut curves indicate that there is negligible acceleration of high— 7 material. Also, the 
functional form of the decrease in < ^v x > is similar for all runs in which 9 = 4° initially. 

Figure 8. The location of the edge of the jet, determined by first finding the location 
along the direction of propagation of the maximum in either 7 ( top curve), fry ( middle 
curve), or p ( bottom curve), then stepping laterally to the point at which the background 
value has increased by 10% of its original value. (The precise value of this constant does 
not affect the results.) The initial maximum radius of the blob is 0.05. After x ~ 1.8 the 
two lower curves coincide. 

Figure 9. The evolution of the weighted mean of the x— component of momentum and 
a measure of the spreading angle <9 >, taken to be the ratio of tangential to axial speeds, 
in the observer’s frame. Shown are (i) < jv x > (solid line) and the reciprocal of < 9 >, 
measured in radians (dashed lines), where the lower limiting 7 value used in the averaging 
for <9 > is taken to be either 4, 5, or 6 ( top panel), (ii) the <9 > values whose reciprocal 
values are indicated in the first panel ( middle panel), and (iii) the number of cells entering 
into the averages for the three 7 cut values in the first two panels ( bottom panel). At late 
times the deceleration makes the < 9 > calculations problematic because the number of 
high 7 cells drops. This progression in the loss of high 7 cells can be seen at t > 4 in the 
bottom panel. 

Figure 10. Light curves constructed by summing the quantity p t ( 1 + f3i cos (f>)/ (l — 

A cos <f>) 2 over the grid. We take this global sum as being a measure of the luminosity seen 
by an observer looking down the jet. The light curves are built up during the course of the 
run by summing the emission from 100 slices moving toward the observer at c that straddle 
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